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We applied a multicanonical algorithm (entropic sampling) to a two-dimensional and a three- 
dimensional Lennard-Jones system with quasicrystalline and glassy ground states. Focusing on the 
ability of the algorithm to locate low lying energy states, we compared the results of the multicanon- 
ical simulations with standard Monte Carlo simulated annealing and molecular dynamics methods. 
We find slight benefits to using entropic sampling in small systems (less than 80 particles), which 
disappear with larger systems. This is disappointing as the multicanonical methods are designed 
to surmount energy barriers to relaxation. We analyze this failure theoretically, and show (1) the 
multicanonical method is reduced in the thermodynamic limit (large systems) to an effective Monte 
Carlo simulated annealing with a random temperature vs. time, and (2) the multicanonical method 
gets trapped by unphysical entropy barriers in the same metastable states whose energy barriers 
trap the traditional quenches. The performance of Monte Carlo and molecular dynamics quenches 
were remarkably similar. 



I. INTRODUCTION 

In the past decade there has been an outpouring 
of interest in accelerating statistical mechanics simula- 
tions. This started with the work of Swendsen and col- 
laborators: Swendsen and Wang introduced a cluster- 
flip method for accelerating non-disordered spin systems 
and Widom, Strandburg, and Swendsen introduced 
a cluster-flip for finding quasicrystalline ground-states in 
a two-dimensional atomic simulation ^ . These methods 
all gain a major speedup by introducing mostly non-local 
update rules, and often prove capable of bypassing criti- 
cal slowing down problems jj]. 

During the same period, different accelerating ap- 
proaches were introduced, which are based on efficient 
schemes to analyze data from traditional Monte Carlo 
simulations and are frequently called "histogram 

methods". These methods have enlarged the applica- 
bility of various kinds of critical phenomena simulations, 
although they are not necessarily designed to bypass crit- 
ical slowing down problems as efficiently as e.g. clus- 
ter algorithms. Nevertheless, substantial progress can 
be achieved combining histogram and cluster-flip algo- 
rithms (see e.g. reference 0). 

More recently so-called "reweighting techniques" have 
been introduced, which are based on an early approach 
by G. M. Torrie and J. P. Valleau They proposed a 
method to enlarge the sampling range of a Monte Carlo 
algorithm by using nonphysical weighting functions. The 
general idea in the newer approaches is to change the rel- 
ative weights of different configurations to sample equally 
in all ranges of energy rather than focusing on a narrow 
temperature range. The most frequently used reweight- 
ing method is multicanonical sampling [p|-^^ , which rep- 
resents the most general method as other reweighting 
methods, e.g. entropic samp ling p^ , can be directly 
mapped onto this approach [|5|. In systems with a 



strongly double peaked probability distribution of mag- 
netization or energy states (a situation often found in 
systems exhibiting a first-order phase transition), the 
multicanonical approach has been proven to be a pow- 
erful tool. Simple reweighting schemes allow to over- 
come the "supercritical slowing down" known from 
canonical Monte Carlo simulations at a fixed tempera- 
ture. E.g. in non-disordered spin systems with a field- 
driven first-order phase transition {e.g. the Ising model) 
or a temperature-driven first-order phase transition {e.g 
the q-state Potts model) the supercritical slowing down 
of canonical Monte Carlo is due to the low Boltzmann 
weight of the domain- wall states. Multicanonical sam- 
pling approaches the problem with introducing a weight- 
function, which weights all magnetization states (Ising 
model) or energy states (Potts model) equally, and there- 
fore ensures that domain- wall states are sampled with the 
same likelihood as all other accessible states. The canon- 
ical distribution function at a fixed temperature, which 
contains all the thermodynamic information, can be re- 
constructed. Usually multicanonical sampling uses local 
update schemes along the lines of the Metropolis algo- 
rithm fl^ ; variations using cluster-flip or other methods 
are feasible and have been proven useful (for a review 
consult ||ll|,|l^ and references therein). 

Of course, acceleration methods are most crucial for 
glassy systems, which otherwise can be inaccessible to nu- 
merical simulations!^. Whether one believes that glasses 
are sluggish because of large energy barriers to relaxation 
(rates ~ e^^"^), or believe that the free energy barriers 
are due to tortuous entropically difficult routes between 



*Even the experiments fall out of equilibrium. Think of an 
experiment as 10^^ parallel atomistic processors with picosec- 
ond clock times! 
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the metastable configurations, a clever algorithm could 
in principle jump directly between the glassy states. In- 
stead of relative rates which grow as a power of T — Tc, 
acceleration could gain us exponential speedups. 

Acceleration methods have been extensively applied to 
disordered spin models. These studies have been less 
focused on understanding the performance of the algo- 
rithms, because the physics of the systems is less thor- 
oughly understood (there has been more to mine from 
the results, and a less firm foundation on which to do 
algorithmic analysis). The multicanonical methods have 
been applied to spin glasses in two and three dimensions 
to calculate the zero-temperature entropy, ground state 
energies, distribution of overlaps etc. (see refs. [|T7|-p5t). 
The authors succeed in evaluating these properties with 
remarkable accuracy; nevertheless, whether the replica 
theory or the droplet scaling ansatz ||2^ is the more 
appropriate picture in describing the ground state prop- 
erties of glasses could not be resolved. The performance 
of multicanonical sampling for glassy systems is clearly 
worse than for system with a less rugged landscape. We 
believe this failure is systematic and can't be avoided 
within the framework of multicanonical sampling. 

The authors of reference jlTj argued multicanonical 
methods should be superior to simulated annealing (grad- 
ual cooling) , although a direct comparison was made 
only to canonical sampling (quenches to a fixed low tem- 
perature). The authors of ref. 29 1 applied multicanonical 
sampling to the Traveling Salesman problem and claim 
to acchieve a dramatic improvement over the traditional 
Monte Carlo simulated annealing approach. Newman 
pO| has used both cluster methods 0| and entropic sam- 
pling to study the random-field Ising model. He finds 
dramatic speedups from both methods, often reaching 
equilibrium in a few passes through the lattice. New- 
man has focused on small systems (mostly 2A^, with a 
few runs for systems up to 64'^), and simultaneously used 
histogram methods to measures critical exponents and 
phase boundaries for a range of disorders and tempera- 
tures. He confirms results of the related "simulated tem- 
pering" approach, invented by E. Marinari and G. Parisi 
pl| , ^ . Simulated tempering proved very useful in spin 
glass simulations and is similar in spirit to the mul- 
ticanonical approach p^ . 

Acceleration methods have been little used in con- 
tinuum atomic simulations, perhaps because of the 
widespread reliance on molecular dynamics methods. 
Straightforward, direct molecular dynamics simulations 
of the equations of motion do not converge to an equi- 
librium state much faster than the Monte Carlo simu- 
lated annealing methods, but they are also not notice- 
ably worse |35| , and they have a direct physical interpre- 
tation. Shumway |Q studied a one-dimensional atomic 
system in an incommensurate sinusoidal potential, and 
developed an evolutionary algorithm which generated op- 
timal cluster moves as the system was quenched to lower 
temperatures; later attempts to generalize these ideas to 
higher dimensions have so far not been successful [p7[. 



The authors of reference used multicanonical sam- 
pling and Monte Carlo simulated annealing to study the 
folding of the peptide Met-enkephalin; the multicanon- 
ical method found the ground state more consistently 
using the same amount of computer time. This result 
underlines the general belief that simulations in 

the multicanonical ensemble are in many ways superior 
to traditional simulated annealing. 

In this paper we apply multicanonical sampling in the 
particular form of entropic sampling to two-component 
Lennard-Jones systems, and compare the performance 
with traditional simulated annealing and straightforward 
molecular dynamics in finding low energy configurations. 
We search for low energy states of a three dimensional 
Lennard-Jones glass, one of the prototype glassy systems 
, and use the set of parameters recently introduced 
by W. Kob and H. C. Andersen [|5|-|7|- addition to 
that, we apply entropic sampling and simulated anneal- 
ing to a two-dimensional Lennard-Jones systems with 
quasicrystalline ground states, using the parameters of 
reference Q . We find that entropic sampling brings little 
benefit for the study of either. We argue that this is likely 
a general effect, applicable to all simulation methods ap- 
plied to glassy systems in the thermodynamic limit. 



II. INTRODUCTION TO THE METHODS: 
MULTICANONICAL SAMPLING AND 
ENTROPIC SAMPLING 

The standard way of implementing a Monte Carlo al- 
gorithm is using importance sampling. The idea behind 
this approach is simple. Rather than weighting each sam- 
ple in phase space equally, each state is weighted with 
a sample probability distribution T(x), where x denotes 
the sampled configuration of the system. To estimate the 
thermal average of an observable A, one calculates: 



j:^Aix)eM-PHix)]T-\x) 
j:^eM-PH{x)]T~\x) 



(1) 



where H is the Hamiltonian of the system (so H{x) is the 
energy E for the state x) and /3 = l/ksT. Choosing r(a;) 
non-uniformly ensures that states with important contri- 
butions to the partition sum are preferentially sampled, 
and therefore the number of states need to be sampled 
to provide a reasonable estimate of A is significantly re- 
duced. 

In standard Monte Carlo methods, i.e. canonical 
Monte Carlo or simulated annealing, the weighting distri- 
bution is the Boltzmann distribution F = exp[—PH{x)]. 
This has the advantage of a direct physical interpreta- 
tion: the computer is doing the same thermal average as 
an equilibrium system at temperature 1/ (fcs/S). It has an 
important disadvantage that configurations and events 
which are rare in the physical system are also rare in the 
simulation. In particular, if the system has a "rugged 
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energy landscape" , with large free energy barriers B sep- 
arating physically important metastable states, the sys- 
tem will cross between these states with the same slow 
Arrhenius rate z/exp(— i?/T) that is found experimen- 
tally. 

The idea of multicanonical sampling is to circumvent 
this problem by choosing r(a;) so that the distribution of 
states P{H{x)) ^ il{H{x)) x r(a;) is approximately 
flat in energy (or some other variable, like magnetization 
1^). In principle, we want to choose r{E) = l/il.{E) = 
exp[— S'(i?)], where Q{E) is the density of states at energy 
E and S{E) is the entropy. Of course, we don't begin the 
simulation knowing the entropy as a function of energy! 

In our work we use the entropic sampling algorithm 
101, which is a numerical and mathematical equivalent 
variant of the multicanonical approach | [l5|] . The only dif- 
ference between entropic and multicanonical sampling is 
the way by which one generates estimates J{E) of S{E). 
The entropic sampling algorithm uses a quite straight- 
forward recursive updating method: 

1: Initialize to zero an array H{E), which will keep 
track of the energy of the visited states. 

2: Sample states according to the current Ji{E), and 
add the energy of each sampled state to the his- 
togram H, for a reasonably long time. 

3: Set the new J^+i according to the following rule: 



J^+l{E) = 



ME)+\og{HiE)), 
ME), 



im{E) 

iiH\E) = 



(2) 



The multicanonical sampling update scheme differs 
from entropic sampling in the treatment of the histogram 
bins with few entries (for an analysis of various schemes 
see references and [^). The original approach in- 
troduces a constant slope for J(E) below a cutoff energy, 
corresponding to a small constant temperature. These 
extra parameters are annoying | |39[ | in the implementa- 
tion; however, they do tend to keep the system from being 
trapped in energy regions which have not hitherto been 
sampled frequently. As we will argue, in glassy systems 
both algorithms will tend to get trapped in low energy 
metastable states even when the statistics are fine[]. In 
this paper we use the simpler entropic sampling method 
of equation (2). 



tin any case our simulations spend around half the time at 
high energies, so any algorithmic improvements can bring at 
best a factor of two in computer time. 



III. THEORETICAL EXPECTATIONS FOR 
RELATIVE PERFORMANCE 

What makes people think that multicanonical sam- 
pling should be an improvement over simulated anneal- 
ing or molecular dynamics? We consider three possible 
reasons. 

(1) Perhaps the multicanonical method is better be- 
cause it allows the system to cross energy barriers (as is 
mentioned frequently fljl)? This is indeed an im- 
provement over canonical sampling at a fixed tempera- 
ture; however, a simulated annealing method also runs 
at a variety of temperatures. 

Indeed, the two methods are identical in the ther- 
modynamic limit. The acceptance ratio for a given 
single-atom Monte Carlo move for entropic sampling is 
P{E) = exp[S'(£;) - S{E')\. In a large system with TV 
atoms, the entropy density S{E)/N is a smooth func- 
tion of the energy density E/N] since the energy den- 
sity change for a single-atom move (£" — E) /N is small, 
we may expand S{E) to first order in E' — E. Us- 
ing the relation dS{E) / dE — the acceptance ra- 
tio becomes P{E) = cx:p[-{E' - E)/T]. Thus entropic 
sampling at the energy E has exactly the same ac- 
ceptance ratio as simulated annealing at a temperature 
T{E) = [dS{E)/dE)-\ 

Thus the local behavior — the acceptance ratio for 
Monte Carlo moves from a given state — is virtually the 
same for multicanonical and canonical sampling. The 
differences between the two methods near a given state 
should be similar in magnitude and type to the differ- 
ences between the microcanonical (fixed-energy) simula- 
tions and the canonical (fixed-temperature) simulations: 
differences can be seen for small systems, but disappear 
as the system gets larger. To be explicit, for a large sys- 
tem the final state of an entropic sampling run for which 
the time-dependent energy is E{t) should be statistically 
equivalent to a simulated annealing run with randomly 
fluctuating temperature T{E(t)). One notes also that the 
quench rate is not tunable for the multicanonical method: 
the "diffusion constant" in energy space depends on the 
atomic step-size and on the number of particles. This is 
potentially a serious handicap, as changing the quench 
rate is the primary tool used in glassy systems to find 
lower energy states. 



^Near first-order transitions, canonical quenches produce 
large changes in the state for small changes in temperature, 
and thus behave quite differently from the multicanonical ap- 
proaches (which by varying the energy explore the interface 
states directly) . This is one of the major applications of multi- 
canonical sampling methods. We expect that the multicanon- 
ical methods will perform for these systems rather similarly to 
microcanonical quenches which conserve the energy; see ref. 
[ ^ . In our glassy simulations, this distinction is presumably 
not important. 
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Thus the power of muhicanonical methods to vary the 
energy to faciUtate barrier crossing is — for large sys- 
tems at least — no different from repeatedly heating and 
cooling the entire system. 

(2) Perhaps the multicanonical method might be pick- 
ing the heating and cooling schedule intelligently, in or- 
der to escape from local minima? Indeed, since the 
effective temperature becomes lower as the energy de- 
creases, an entropic sampling system stuck in a high- 
energy metastable state will have larger thermal excita- 
tions (bigger acceptance of upward moves in energy) than 
one in a low energy state, and will depart faster. This 
is the explanation, we believe, for the substantial success 
of the entropic (multicanonical) sampling method seen in 
the past. 

This preferential escape from high-energy metastable 
states will unfortunately also become unimportant for 
large systems. One can see this most easily by consider- 
ing a local region trapped in a high-energy configuration 
with local energy e', and with a lower energy configura- 
tion e nearby, separated by a barrier b. For a small sys- 
tem, where the local energy difference e' — e is important, 
the effective temperatures in states e' and e will differ, 
but for a large system of size N this temperature differ- 
ence (from the differing acceptance ratios from the two 
states) will vanish as 1/N. There are glassy systems in 
which the energy barriers and energy differences are not 
all local: mean-field spin glasses, for example, have en- 
ergybarriers which grow as powers of the number of spins 
N ||2q| . However, the maximum energy barrier (and pre- 
sumably the maximum energy asymmetry e' — e) scales 
with a power TV" with a strictly less than one (at least in 
finite dimensions), so the change in effective temperature 
N"/N stiU vanishes as ^ oo ||l|l. 

(3) Perhaps the multicanonical method is exploring dif- 
ferent energies more effectively than an externally chosen 
cooling schedule for simulated annealing? For example, 
the multicanonical method is guaranteed to converge to 
an equilibrium density of states at each energy. The same 
is true for a simulated annealing run at infinitely slow 
cooling, but is not true for repeated coolings at a fixed 
rate, which would be expected to generate metastable 
states repeatedly. On the one hand, the theorems sug- 
gest that the ground state should be occupied as often 
as any other energy; on the other hand it is hard to see 
how a multicanonical quench to low energies can bypass 
the metastable states that trap simulated annealing runs 
of comparable computer time. 

To address this issue, let us consider the characteris- 
tics of the random walk E{t) that the system performs 
in energy space as a function of time within the multi- 
canonical approach. For some systems, such as the Ising 
model, multicanonical sampling does indeed produce a 
roughly unbiased random walk (if one starts from a good 
estimate of the density of states). As the system becomes 
larger, the energy range scales with N and the step-size 
of the energy stays fixed, so the time-scale for diffus- 
ing from high energies to near the ground-state scales 



as N"^ (distance scales with the square root of time). 
This behavior is confirmed in simulations [^,0 study- 
ing, e.g. the first-order phase transition below Tc chang- 
ing the Ising model from pointing up to down, where 
traditional canonical methods would suffer from the sur- 
face tension barrier aL'^^^ — crA^t''"!)/'' and so the time 
scales as expaN'-'^"^^^'^. Bypassing this "supercrit- 
ical slowing-down" [p^ is an important application for 
multicanonical methods. 

It is known numerically that this simple argument 
breaks down in simulations of spin glasses . The typ- 
ical time to cover the energy range (called the ergodicity 
or tunneling time in the literature) for spin glasses scales 
as N'^ ^ or perhaps Q instead of iV^. Why should 
the random walk argument not work for glasses? 




Temperature or Energy 



FIG. 1. The glassy metastable states are often thought to 
form a tree-like structure. We anticipate that the inaccessible 
metastable states will form a barrier to leaving the ground 
state within the multicanonical approach, because the total 
density of states grows much faster than the accessible density 
of states. 

The answer to why the random walk argument breaks 
down, we assert, can be found in the trapping of entropic 
sampling due to the same metastable states explored in 
thermal coolings. Indeed, it is the metastable states that 
the system is not able to explore that trap the entropic 
sampling algorithm. The states of glassy systems are 
often described in a caricature tree-like structure (side- 
ways in figure 1). The horizontal axis of the tree can be 
thought of as either energy or temperature: the branches 
represent mutually inaccessible ergodic components. For 
the Ising model, there are two major ergodic components 
(corresponding to the two directions for the magnetiza- 
tion) and a few domain-wall states. For glasses, the er- 
godic components are sometimes thought of as regions 
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of configuration space separated by infinite free energy 
barriers (as in the mean- field spin glass models p6| ), and 
sometimes as regions separated by energy barriers which 
are too large to cross in the time-scale of the experiment 
or simulation. 

The key point is that the accessible density of states 
for a glass can be very different from the total density 
of states. In figure I, we note that the ergodic compo- 
nent containing the ground state has a density of states 
which differs from the density of states for the system 
as a whole, starting at the energy of the first accessible 
metastable states. The number S of these inaccessible 
metastable states is related to the density of tunneling 
states in configurational glasses [^-p^, and is thought 
to increase exponentially with the size of the system (M 
independent two-state systems with uncrossable barriers 
generate E — 2*^ states). In spin glasses, the number 
of components separated by infinite free energy barriers 
(ones which diverge as ^ oo) diverges with a power 
ofiV 




.Ground State Well K ~ 
Temperature or Energy 



FIG. 2. Grossly oversimplified picture of of glassy states, 
assuming that all components merge at the glass transition 
at energy Eg. The parameter ^ denotes the fraction of er- 
godic components at each energy or temperature. At the 
glass transition, the fraction of accessible ergodic components 
is of order thus the time to find the ground state ergodic 
component (or leave it once found) scales with E as S — * oo. 

The multicanonical sampling method is guaranteed to 
sample the ground-state energy just as much as any other 
energy. It is easy to see, however, that once the system 
is in the ground state, it will stay there for a long time! 
If the density of states within the ground-state ergodic 
component is Qg{E) and the density of states in the en- 
tire system is ^{E), then the acceptance ratio for a mul- 
ticanonical sampling move from E to E' is Q.{E) / Vi{E' ) , 



while the probability of a random move raising the en- 
ergy is Vlg{E')/Clg{E). Hence the likelihood of sampling 
high energy states E within the ground-state component 
will fall as Q.g{E) /Vl{E). Consider the very crude model 
where all ergodic components are similar and stay com- 
pletely separate until the glass transition energy Eg (the 
energy at the glass transition temperature), at which 
point they merge (figure 2). For this model, escaping 
from the ground state component will take a time which 
scales as the total number of ergodic components, and 
hence diverges as N — > oo. Since multicanonical sam- 
pling spends the same amount of time in each energy 
range, the time between independent visits to the true 
ground state will scale as the time needed to escape from 
the ground state ergodic component. 

So, we begin our exploration with the expectation that 
multicanonical sampling should be useful for small sys- 
tems, but will not provide significant advantages for large 
system sizes. 

IV. IMPLEMENTATION 

We argued in the previous section that entropic sam- 
pling will not be a fundamental improvement over re- 
peated coolings using simulated annealing, at least in 
large systems. On the other hand, entropic sampling 
and other multicanonical methods have been reported to 
lead to substantial gains in equilibration times for small 
glassy spin systems |L7|-|25| . We are interested in simula- 
tions of structural glasses: collections of atoms which typ- 
ically form metastable, glassy configurations when slowly 
cooled. In this section, we give a detailed description of 
our implementation of entropic sampling, simulated an- 
nealing, and molecular dynamics. In order to ensure a 
fair comparison, we have as far as possible taken cooling 
schedules and time and spatial step sizes from standard 
references in the literature. 

In our three dimensional simulations, we applied the 
three algorithms to a binary mixture of large (L) and 
small (S) particles with the same mass, interacting 
via the Lennard- Jones potential of the form Vap{r) = 
4:eai3[{craf:i/ry^ — (cTap/r)^]- The values of e and a were 
chosen as follows: eLL — ^-O^o-ll = 1-0, cls — 1.5, (Tls — 
0.8, ess = 0.5, ass = 0.88. All results are given in re- 
duced units, where aLL was used as the length unit and 
Ell as the energy unit. The systems were kept at a fixed 
density (p w f-2), periodic boundary conditions have 
been applied and the potential has been truncated appro- 
priately according to the minimum image rule [^5| , and 
shifted to zero at the respective cutoff. The minimum 
image rule prevents a particle from using the periodic 
boundary conditions to see more than one copy of its 
neighbors: to use the conventional cutoff at r = 2.5(7 
would demand at least 160 particles. The choice of 
parameters follows recent simulations of Lennard Jones 
glasses [p5[-^ ; this choice suppresses recrystallization of 
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the system on molecular dynamics time scales. This po- 
tential together with this set of parameters mimics the 
potential for Ni8oP2o- We looked at 5 different system 
sizes (A^=20, 40, 60, 80, 100). For each N we gener- 
ated 30 low energy configurations. The initial configu- 
rations were random in the case of simulated annealing 
and entropic sampling, and high temperature equilibrium 
configurations in the molecular dynamics case. To com- 
pare the three methods, we defined one run length to be 
10^ sweeps through the system for the two Monte Carlo 
methods. The molecular dynamics runs are quenched at 
a rate which consumes the same CPU time as used by 
the Monte Carlo sampling. 

In our two dimensional simulations, we did not apply 
molecular dynamics (our hard-wall boundary conditions 
made it inconvenient). We again used a binary Lennard 
Jones system, introduced by Widom, Strandburg, and 
Swendsen with a slightly unconventional form for the 
potential: V^pir) = e„;3[(CT„^/r)i2 - 2{aap/rf]. The 
Lennard Jones parameters are chosen to favor configura- 
tions of decagonal order (e^L = 0.5, (7^^ = 1.176, e^g = 
1.0, a^s = l-Oi ess = 0.5, css = 0.618), and the sys- 
tem is known to have a quasicrystalline ground state. 
The particles are initially randomly distributed in a large 
cylindric box with infinitely high walls. The potential 
was truncated at r cutoff = ^.baap and shifted to zero 
at this point. All results are in reduced units with e^s 
and <7ls fundamental units. Here 4 different system 
sizes (N=31, 66, 101, 160) were used, where the num- 
ber of different particles were chosen to keep the ratio 
fixed close to the value of 1.06 large atoms per small 
atom. The authors of reference found that this ra- 
tio led to defect-free ground states. This system pro- 
vides an excellent testing ground for entropic sampling 
for various reasons. The ground state is known to be 
quasicrystalline, a state with a strong bond orientational 
order without a long-range order periodicity. Defective 
configurations are easily recognized, as the typical defects 
consists of triangles of like particles. There are plenty of 
metastable states with high energy barriers, as it takes 
rearrangement of a large number of particles to disen- 
tangle the triangular defects. The authors of reference 
1^ have shown that simulated annealing fails to locate 
ground states, and always gets trapped in a long lived 
metastables states, a problem they circumvented using 
three-particle cluster-flips. 

The final minimum energy configurations for runs of all 
three methods were optimized by starting from the lowest 
energy configurations found, and quenching down to T = 
using a conjugate gradient method. The resulting T — 
energies are compared in the following section. 

The entropic sampling method was implemented gen- 
erally as described in section 2. The Metropolis algo- 
rithm [|l6| was used for local updates. We developed an 
initial estimate J{E) for the entropy S{E) with a long 
run (approximately 10^ sweeps) starting from a flat dis- 
tribution: this initial estimate was used as the starting 
distribution for the subsequent runs. In three dimensions 



we redid this initialization for each system size; in two di- 
mensions we initialized in this way for the 101 particle 
system and used finite-size scaling [|l^ from this distri- 
bution for initialization at other system sizes. Finite-size 
scaling does not appear to work for glassy systems with 
quenched disorder . We did not include this initial 

computer time in the comparisons: thus we err on the 
side of entropic sampling. Notice that in continuous sys- 
tems it is not obvious how to set the optimal bin size for 
the histogram (unlike in spin systems, where the small- 
est energy step determines the bin size). We therefore 
tested several bin sizes. For the two-dimensional systems 
a fixed bin size of 0.001 has been used, and in the three- 
dimensional systems, we found it useful to set the bin size 
to 0.01/A^ energy units. Our investigations suggest that 
larger bin sizes can introduce artificial barriers in the 
low energy range, and smaller bins lead to more noise. 
To set a context, the typical successful energy step in 
a 20 particle simulation in three dimensions varied from 
around one at high temperatures to around 0.01 near the 
ground state. We explored energy-dependent step sizes 
for the single-atom moves, but they did not improve per- 
formance. Entropic sampling demands an upper cutoff 
for the energy: we use zero for the upper limit for both 
two and three dimensions. 

Simulated annealing uses locally the Metropolis update 
scheme with the Boltzmann factor as the sample prob- 
ability distribution. The cooling schedule implemented 
here is similar to the one used in references P,|38[| , where 
the temperature is repeatedly lowered by a small factor 
and then annealed. In our runs, we choose fifty annealing 
steps of 20,000 sweeps each, with an initial temperature 
of one and a final temperature of 0.05; each tempera- 
ture is thus cooled down by a factor of 0.942. The initial 
configurations were set at random. 

The molecular dynamics routine used the velocity 
form of the Verlet algorithm |5^. The unit of time 
is given by {maj^ / 4:8e l lY , where m is the mass 
of the particle: the Verlet time step in these units is 
6t = 0.01. The system was coupled to a heat bath and 
the temperature was reduced linearly in time according 
to Tbath = Tltart " 7md X t. Note that the cooHng 
here is linear in time , as is traditional in molecular dy- 
namics of Lennard- Jones glasses . The cooling rate 
7md — 1-0 10"'' was chosen so that the MD runs consume 
an amount of computer time similar to that of the Monte 
Carlo algorithms. This cooling rate is in the middle of 
the range explored in recent simulations, although our 
system sizes are much smaller psf . The initial configura- 
tions were equilibrated at a temperature Tgtart = 1.0 (at a 
small cost of computer time which we did not factor into 
the comparisons), and cooled to the final temperature 
Tfinai = 0.05, yielding approximately 2.0 x 10^ molecular 
dynamics steps. 



6 



V. RESULTS 

In this section we will first compare the performance 
of the three methods in locating low energy states of the 
two and three-dimensional Lennard- Jones systems. The 
performance of the three methods is remarkably similar. 
Second, we will compare low energy configurations of the 
two-dimensional system to show that the the algorithms 
get trapped in similar metastable states. Third, we will 
quantitatively analyze the trapping of the entropic sam- 
pling algorithm in a metastable state. 

We present the T — energies of the lowest energy 
configurations for the three-dimensional Lennard-Jones 
systems in Table 1. For = 20 particles each algorithm 
is able to locate the same lowest energy state, presumably 
the ground state. For iV = 40 and = 60 particles the 
lowest energy state is found by entropic sampling. The 
gain in energy AE^q over simulated annealing is around 
0.03, and the gain is 0.02 over molecular dynamics. For 
= 80 and N — 100 particles the lowest energies are 
found by molecular dynamics, and the gain over entropic 
sampling is AE^o ~ 0.05 and AEiqq ^ 0.02. 



TABLE 1: T = energy per particle for the lowest energy con- 
figuration found with entropic sampling, simulated annealing and 
molecular dynamics. 



N 


Entropic 


Simulated 


Molecular 




Sampling 


Annealing 


Dynamics 


20 


-0.89 


-0.89 


-0.89 


40 


-4.18 


-4.15 


-4.16 


60 


-5.38 


-5.36 


-5.37 


80 


-6.52 


-6.56 


-6.57 


100 


-6.85 


-6.86 


-6.87 



There are three things to notice about this table. First, 
the dramatic energy difference with increasing system 
size is due to the change in the cutoff in the potential 
given by the minimum image rule. Using the twenty 
particle cutoff in the larger system sizes, we found en- 
ergies which hardly varied with system size. Second, the 
fact that these energies differ in the third decimal place 
does not mean that the differences are negligible. In ref. 
p8| , H the dependence of the final energy on the cooling 
rate for exactly this system was studied using molecular 
dynamics: to gain an energy of 0.03 starting from the 
cooling rate we are using, they had to decrease the cool- 
ing rate by a factor of ten. Third, as we argued in section 
three, any gains given by entropic sampling disappear as 
the system size grows. 

In Table 2 we list the mean and the standard deviation 
of the energies from the thirty runs at each system size 
with each algorithm. It has been found in the literature 
that the fiuctuations for simulated annealing are much 
larger than for entropic sampling We find this to 

be true both for simulated annealing and for molecular 
dynamics. Indeed, the average performance of entropic 
sampling remains comparable to that of the other two 



methods, even at the larger system sizes (where the ex- 
tremal performance was worse). 



TABLE 2: Mean energy per particle and the standard deviation 
evaluated using all low energy configurations found by entropic 
sampling, simulated annealing and molecular dynamics. 



N 


Entropic 


Simulated 


Molecular 




Sampling 


Annealing 


Dynamics 


20 


-0.84 ± 0.03 


-0.79 ± 0.06 


-0.81 ± 0.04 


40 


-4.10 ± 0.03 


-4.07 ± 0.05 


-4.08 ± 0.04 


60 


-5.34 ± 0.01 


-5.33 ± 0.03 


-5.33 ± 0.03 


80 


-6.50 ± 0.01 


-6.51 ± 0.02 


-6.51 ± 0.03 


100 


-6.83 ± 0.01 


-6.84 ± 0.02 


-6.82 ± 0.02 



The Holy Grail of this field is to accelerate three- 
dimensional glass simulations, bypassing barriers to re- 
laxation. Perhaps this is too high a standard — nobody 
has such an algorithm. We now will apply entropic sam- 
pling to a two-dimensional Lennard- Jones system, where 
an effective cluster-flip acceleration method has been de- 
veloped Hj. In Tables 3 and 4 we show the extremal 
and the mean T = energies for a variety of system 
sizes. Again, entropic sampling is slightly better for the 
smaller systems, but the advantage disappears for the 
largest system. 



TABLE 3; Energy per particle for the lowest energy configuration 
found with entropic sampling and simulated annealing. 



N 


Entropic Sampling 


Simulated Annealing 


31 


-2.1 


-2.1 


66 


-2.27 


-2.22 


101 


-2.34 


-2.33 


160 


-2.37 


-2.38 



TABLE 4; Mean energy per particle and the standard deviation 
evaluated using all low energy configurations found by entropic 
sampling and simulated annealing. 



N 


Entropic Sampling 


Simulated Annealing 


31 


-1.97 ± 0.1 


-1.87 ± 0.17 


66 


-2.24 ± 0.04 


-2.12 ± 0.07 


101 


-2.29 ± 0.07 


-2.23 ± 0.06 


160 


-2.34 ± 0.05 


-2.34 ± 0.02 



It is remarkable how similarly the three different meth- 
ods perform. Although it seems to be known that molec- 
ular dynamics and simulated annealing are comparable 
psf , we are not aware of any reference providing a di- 
rect comparison. Of course, comparisons of efficiency 
are highly implementation dependent. The two Monte 
Carlo methods could benefit from a temperature depen- 
dent step-size (although we did experiment with it with- 
out finding any substantial improvement). One could 
refine the cooling schedule for the two traditional meth- 
ods. One could introduce a temperature cutoff (like in 
the original multicanonical approach) or use variable bin- 
sizes to improve the entropic sampling method. Again, 
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our experiments with bin-size and cutoff were not en- 
couraging. Our main conclusion is that the choice of 
methods is a matter of taste. In particular we are encour- 
aged by the fact that Monte Carlo methods are competi- 
tive, especially as the adapt easily to cluster acceleration 
methods. 

All three methods suffer from the large number of 
metastable states prevalent in the configuration space of 
the two and three-dimensional systems, and thus are not 
capable of locating ground states. To show that they find 
similar metastable states, we plot in figures 3 and 4 the 
lowest energy configuration found by entropic sampling 
and by simulated annealing. 




FIG. 3. Lowest energy configuration generated with en- 
tropic sampling. The grey particles forming a triangle repre- 
sent a defect. 




FIG. 4. Lowest energy configuration generated with sim- 
ulated annealing. The grey particles forming a triangle repre- 
sent a defect. 

These configurations are typical representatives of 
metastable states for the two-dimensional system. The 



defects are clusters of three large particles, which are 
shown in grey in figures 3 and 4. 

Why is entropic sampling not bypassing the free- 
energy barriers to relaxation? Wc finish this section with 
a vivid illustration of how the entropic sampling algo- 
rithm gets trapped in a metastable state. In the bottom 
half of figure 5 we plot the energy as a function of time: 
at very short times it performs a random walk in energy 
space as advertised, but it rapidly gets trapped in a low 
energy metastable state. 

The simulation shown in figure 5 is the same as the 
runs for N = 100 particles tabulated in Tables 1 and 
2 except for two important differences. (1) The runs of 
Table 1 and 2 ran for 10^ sweeps, here we ran for 10^ 
sweeps. (2) The entropy estimate J{E) for the runs in 
Tables 1 and 2 was dynamically updated every 10^ sweeps 
using the recursive updating scheme equation (2). Here 
we calculated a best estimate J{E) from the thirty runs 
in the Tables and used this functi on as a fixed entropy 
estimate. The best estimate J{E) (comprising informa- 
tion from 40 x 10^ sweeps) is a sufhciently smooth func- 
tion that we don't expect (or observe) the system to be 
trapped in some artificial well resulting from statistical 
fluctuations in J{E). 



20.0 




-6.8 -5.8 -4.8 -3.8 -2.8 -1.8 -0.8 
Energy per particle 



FIG. 5. An entropic sampling run, showing trapping into 
a metastable state. In this run, we used a fixed J{E) {i.e. no 
dynamical updating) gleaned from several previous runs of 
the same system. Energy vs. time is shown on the bottom 
panel. Only the first quarter of the time series is shown. At 
short times, we observe a random walk in energy. At times 
below 2 X 10^ sweeps the system randomly walks through ener- 
gies above the glass transition. The system then gets trapped 
in two successively lower metastable free-energy wells (see 
figure 6). The Histogram log{H{E)) = AJ{E) of visited 
states measured during the simulation is shown in the up- 
per panel. Note the large peak eissociated with the trapping. 
The second and third peak are above the glass transition and 
correspond to more transient states. 
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The top half of figure 5 shows the logarithm of the 
histogram H{E) tabulating the visited states as a func- 
tion of energy. This function is important as it is used 
in the recursive updating scheme \og{H{E)) = AJ{E) = 
J {E) update - -'^(£^)cstimatc (sse equation (2)). 

Figure 6 shows that the system is trapped in a single 
harmonic metastable state. The upper panel shows an 
expanded view of the first peak in AJ{E). For times 
after 1.3 x 10^ the system exclusively samples in a single 
well: repeated quenches yield the same minimum energy 
El = —6.8241. This is a metastable state: as seen in 
Table 1 the true ground state has an energy < —6.87. 



20.0 




-6.80 -6.75 -6.70 

Energy per particle 



FIG. 6. Expanded view of figure 5 at low energies. En- 
ergy vs. time is shown on the bottom panel. From 2 x 10^ 
to 1.3 X 10® sweeps, the system seems trapped in a num- 
ber of states: repeated quenches yield metastable energies 
clustering around E2 = —6.82 ± 0.002. The system then 
falls into a slightly lower state, with energy Ei — —6.8241; 
repeated quenches show that it stays in this single well 
for the remainder of the simulation. The equilibration 
within this lowest well is excellent: the acceptance ratio 
is near 50% and the system exhibits an efficient random 
walk in energy within the single well. The Histogram 
\og{H{E)) = AJ{E) of visited states is shown in the upper 
panel as a dark line. The light line is the theoretical prediction 
assuming a single, harmonic well at energy Ei = —6.8241: 
AJharmonic(£;) = (3iV/2 - 1) log(S - El) - J(E) + C. The 
bumps in the theoretical curve are due to the irregularities 
in our initial estimate J{E); the identical-looking bumps in 
the measured data reflect the efi'ects of J{E) in weighting the 
histogram. Our theoretical prediction describes the measured 
data well for the energy range explored during the last por- 
tions of the simulation. The range above E/N ~ —6.78 is 
underestimated as in this region the data is composed largely 
from states corresponding to the minima clustered around £2. 

In the harmonic approximation we can analytically cal- 
culate the density of states and compare the contribution 



from the metastable state directly to the measured data. 
The harmonic density of states has the form 

3N 1 

K j ' 

where K involves the geometric mean of the phonon fre- 
quencies and can be thought of as a typical spring con- 
stant. In the entropic sampling algorithm the probability 
of sampling a state in this harmonic well is given by the 
ratio of the density of states in the s ingle w ell divided by 
the estimated density of states exp(J(£')). This proba- 
bility is compared directly to the histogram of sampled 
states in the upper half of figure 6. The system is trapped 
in a single harmonic well. 

By running for shorter times and by dynamically up- 
dating the entropy estimate during each run, we have 
substantially mitigated the trapping problem of entropic 
sampling shown in figures 5 and 6. Imagine the entropy 
estimate J{E) after updating it by adding \og{H{E)). 
The acceptance ratio to leave the region given by 1/ J{E) 
will dramatically increase and thus the trapping will be 
bypassed, as indeed we observed in practice. Lingering 
near a state increases the estimate entropy in that region 
and eventually push the system out. But dynamical up- 
dating should not be an essential ingredient of the algo- 
rithm, which is formulated presuming an a priori knowl- 
edge of the entropy as a function of energy. Formally 
dynamical updating violates the Markovian character of 
the algorithm and convergence to the equilibrium state 
is no longer guaranteed. In practical terms it is very dis- 
tressing that the algorithm needs to produce a noticeable 
bump in the density of states to escape from a metastable 
state. The comparison against molecular dynamics and 
simulated annealing would be substantially more unfa- 
vorable for long runs without dynamical updating. 

The figures 5 and 6 are a tangible illustration of the 
trapping mechanism depicted in figures 1 a nd 2. T he in- 
accessible metastable states contributing to J{E) form a 
strange type of entropic barrier around the metastable 
state £^1. Leaving Ei via a saddlepoint at E2 > 
Epcak ^ —6.8 is suppressed by roughly the exponential 
of AJ{E2) — AJ(Epcak) shown in figure 6. 

Slow cooling in molecular dynamics or simulated an- 
nealing can lead to trapping in metastable states due to 
large energy barriers. Entropic sampling and the other 
multicanonical methods get trapped in metastable states 
because of large entropic barriers imposed by the algo- 
rithm. In both case the algorithms are sabotaged by the 
large number of low lying metastable states. Entropic 
sampling provides a new insight into this problem but 
doesn't provide a solution. 

VI. CONCLUSIONS 

In this study we applied the multicanonical method en- 
tropic sampling to Lennard-Jones systems. We focused 
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on the ability of the algorithm to find ground states of 
these glassy systems and compared the performance to 
the two traditional glassy simulation methods simulated 
annealing and molecular dynamics. The use of entropic 
sampling didn't reveal any new insights into the ground 
state properties of Lennard-Jones glasses. We explain 
these results on the basis of the following observations. 

First, in the thermodynamic limit multicanonical 
methods are locally equivalent to simulated annealing. 
Furthermore the global dynamics of multicanonical sam- 
pling resembles a random heating and cooling of the sam- 
ple. Thus for large systems simulated annealing and mul- 
ticanonical sampling must have the same properties. In 
principle multicanonical sampling has the advantage of 
providing the density of states, which allows to evalu- 
ate the canonical distribution function. In glasses this 
feature is not necessarily helpful, as the multicanonical 
methods samples phase space as slowly as the annealing 
methods, thus in practice multicanonical sampling will 
not be able to extract any equilibrium expectation val- 
ues better than simulated annealing. 

Second, the large number of inaccessible metastable 
states imposes a bizarre entropy barrier to the multi- 
canonical method. The algorithm simply gets stuck in a 
metastable state, as it might using molecular dynamics 
and simulated annealing. We underlined this point by 
comparing the probability distribution estimated by the 
algorithm inside the metastable state with a theoretical 
expression derived in the harmonic approximation. 

Furthermore our results emphasize the known fact, 
that simulated annealing and molecular dynamics have 
similar performance in glassy systems. As a consequence 
one should acknowledge the importance of averaging 
over many molecular dynamics trajectories especially for 
glassy systems. Averages over an ensemble of trajecto- 
ries are a basic concept in Monte Carlo simulations, the 
striking similarity in performance to molecular dynamics 
simulations is a hint to the importance of similar averages 
in glassy molecular dynamics simulations. 

Finally, the goal of finding a method which gains an 
exponential speed-up of glassy simulations still remains. 
Our study clearly indicates that standard reweighting 
techniques will presumably be no substantial help in 
tackling this problem. The complicated structure of the 
glassy configuration space needs more intelligent algo- 
rithms, which are not only able to bypass energy barriers 
but also to find an efficient path through the rugged en- 
ergy landscape. 
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